#!/usr/bin/env julia
"""
Diagnostic script to dump wind BUFR metadata for parity comparison.

This script reads wind observations using BUFRLIB and dumps detailed metadata
to compare with Fortran GSI fort.207 output.

Usage:
    julia --project=. scripts/julia/diagnostics/dump_wind_bufrlib_metadata.jl [prepbufr_file] [output_csv]

Arguments:
    prepbufr_file: Path to PrepBUFR file (default: tutorial case)
    output_csv: Output CSV path (default: results/wind_bufrlib_metadata.csv)

Output columns:
    index, station_id, longitude, latitude, pressure_Pa, u_ms, v_ms,
    platform_code, report_type, surface_cat, quality_mark, sigma_mult,
    duplicate_factor, anemometer_type, obs_height_m, elevation_m, time_offset_hr
"""

using Dates
using Printf
using Statistics
using CSV
using DataFrames

# Add src to load path
push!(LOAD_PATH, normpath(joinpath(@__DIR__, "..", "..", "..", "src")))

using GSICoreAnalysis.DataIO.WindBUFRReader

function dump_wind_metadata(prepbufr_file::String, output_csv::String;
                           max_obs::Int = 100000,
                           verbose::Bool = true)
    println("=" ^ 80)
    println("Wind BUFR Metadata Diagnostic")
    println("=" ^ 80)
    println("Input file: $prepbufr_file")
    println("Output CSV: $output_csv")
    println()

    # Configure reader
    config = WindBUFRConfig(
        use_bufrlib = true,
        max_obs = max_obs,
        qc_threshold = 15,  # Accept all QC levels for diagnostic
        apply_10m_reduction = false,  # Keep raw values for comparison
        verbose = verbose
    )

    # Read observations
    println("Reading wind observations with BUFRLIB...")
    bundle = read_wind_bufr(prepbufr_file, config)

    println("Read $(bundle.nobs) wind observations")
    println()

    # Summary statistics
    println("Summary Statistics:")
    println("-" ^ 40)

    platform_counts = Dict{Int,Int}()
    qc_counts = Dict{Int,Int}()
    for obs in bundle.observations
        platform_counts[obs.platform_code] = get(platform_counts, obs.platform_code, 0) + 1
        qc_counts[obs.quality_mark] = get(qc_counts, obs.quality_mark, 0) + 1
    end

    println("Platform distribution:")
    for (code, count) in sort(collect(platform_counts), by=x->x[2], rev=true)
        platform_name = get(WindBUFRReader.READ_WIND_PLATFORM_MAP, code, :unknown)
        pct = 100.0 * count / bundle.nobs
        @printf("  %3d (%-12s): %6d obs (%5.1f%%)\n", code, platform_name, count, pct)
    end

    println("\nQuality mark distribution:")
    for (qc, count) in sort(collect(qc_counts))
        pct = 100.0 * count / bundle.nobs
        @printf("  QC %2d: %6d obs (%5.1f%%)\n", qc, count, pct)
    end

    # Wind speed statistics
    speeds = [sqrt(obs.u^2 + obs.v^2) for obs in bundle.observations]
    println("\nWind speed statistics (m/s):")
    @printf("  Mean:   %.2f\n", mean(speeds))
    @printf("  Median: %.2f\n", median(speeds))
    @printf("  Min:    %.2f\n", minimum(speeds))
    @printf("  Max:    %.2f\n", maximum(speeds))
    @printf("  Std:    %.2f\n", std(speeds))

    # Pressure level distribution
    pressures_hpa = [obs.pressure / 100.0 for obs in bundle.observations]
    println("\nPressure level statistics (hPa):")
    @printf("  Mean:   %.1f\n", mean(pressures_hpa))
    @printf("  Median: %.1f\n", median(pressures_hpa))
    @printf("  Min:    %.1f\n", minimum(pressures_hpa))
    @printf("  Max:    %.1f\n", maximum(pressures_hpa))

    # Check for extended metadata availability
    has_dup_factor = count(obs -> obs.duplicate_factor != 1.0, bundle.observations)
    has_anemo_type = count(obs -> obs.anemometer_type != 0, bundle.observations)
    has_height = count(obs -> obs.observation_height > 0.0, bundle.observations)

    println("\nExtended metadata availability:")
    @printf("  Duplicate factors (≠1.0): %6d obs (%5.1f%%)\n",
            has_dup_factor, 100.0 * has_dup_factor / bundle.nobs)
    @printf("  Anemometer types (≠0):    %6d obs (%5.1f%%)\n",
            has_anemo_type, 100.0 * has_anemo_type / bundle.nobs)
    @printf("  Observation heights (>0): %6d obs (%5.1f%%)\n",
            has_height, 100.0 * has_height / bundle.nobs)

    println()

    # Create DataFrame for CSV output
    df = DataFrame(
        index = 1:bundle.nobs,
        station_id = [obs.station_id for obs in bundle.observations],
        longitude = [obs.longitude for obs in bundle.observations],
        latitude = [obs.latitude for obs in bundle.observations],
        pressure_Pa = [obs.pressure for obs in bundle.observations],
        u_ms = [obs.u for obs in bundle.observations],
        v_ms = [obs.v for obs in bundle.observations],
        platform_code = [obs.platform_code for obs in bundle.observations],
        report_type = [obs.report_type for obs in bundle.observations],
        surface_cat = [obs.surface_category for obs in bundle.observations],
        quality_mark = [obs.quality_mark for obs in bundle.observations],
        sigma_mult = [obs.sigma_multiplier for obs in bundle.observations],
        duplicate_factor = [obs.duplicate_factor for obs in bundle.observations],
        anemometer_type = [obs.anemometer_type for obs in bundle.observations],
        obs_height_m = [obs.observation_height for obs in bundle.observations],
        elevation_m = [obs.elevation for obs in bundle.observations],
        time_offset_hr = [obs.time_offset for obs in bundle.observations]
    )

    # Write CSV
    CSV.write(output_csv, df)
    println("Wrote metadata to: $output_csv")
    println()

    # Sample observations for visual inspection
    println("Sample observations (first 10):")
    println("-" ^ 40)
    for i in 1:min(10, bundle.nobs)
        obs = bundle.observations[i]
        platform_name = get(WindBUFRReader.READ_WIND_PLATFORM_MAP, obs.platform_code, :unknown)
        @printf("%3d: %-8s (%.2f, %.2f) P=%6.1f hPa u=%6.2f v=%6.2f platform=%d (%s) QC=%d\n",
                i, obs.station_id, obs.longitude, obs.latitude,
                obs.pressure/100.0, obs.u, obs.v, obs.platform_code, platform_name, obs.quality_mark)
    end

    println()
    println("=" ^ 80)
    println("Diagnostic complete!")
    println("=" ^ 80)

    return bundle, df
end

# Main execution
function main()
    # Default paths
    default_prepbufr = "/home/docker/comgsi/tutorial/comGSIv3.7_EnKFv1.3/run/2018081212/obs/rap.t12z.prepbufr.tm00"
    default_output = "results/wind_bufrlib_metadata.csv"

    # Parse arguments
    prepbufr_file = length(ARGS) >= 1 ? ARGS[1] : default_prepbufr
    output_csv = length(ARGS) >= 2 ? ARGS[2] : default_output

    if !isfile(prepbufr_file)
        println("ERROR: PrepBUFR file not found: $prepbufr_file")
        println()
        println("Usage: julia dump_wind_bufrlib_metadata.jl [prepbufr_file] [output_csv]")
        exit(1)
    end

    # Ensure output directory exists
    mkpath(dirname(output_csv))

    # Run diagnostic
    dump_wind_metadata(prepbufr_file, output_csv)
end

if abspath(PROGRAM_FILE) == @__FILE__
    main()
end
